library(foreign)
library(ergm)
library(xtable)
library(combinat)
library(ROCR)

######### Data Preoparation ###########
# Read in Atop alliance data - 1816-2001
atopdy <- read.csv("./Data/source_data/atop3_0dyeug.csv")

# Read in state system-year data
ssys <- read.csv("./Data/source_data/system2008.1.csv")

# Read in state bound dates
stdates <- read.csv("./Data/source_data/states2008.1.csv")

# Read in CINC data - 1816-2000
cinc <- read.csv("./Data/source_data/CINC.csv")

# Read in Maoz dataset (use only trade seq) - 1816-2000
maoz <- read.dta("./Data/source_data/jcr2006replication.dta")

# Read in contiguity data (one record per contiguity relationship) - all yrs
contig <- read.csv("./Data/source_data/contdir.csv")

# Polity - 1800-2007
polity <- read.csv("./Data/source_data/polity.csv")

# Read in MIDs Data (use maoz)
mid <- read.csv("./Data/source_data/MIDs.csv")

# Read in Adverse Regime Changes Data - 1952-2007
advch <- read.csv("./Data/source_data/Adverse_Chg55.csv")

# Read in major powers
majp <- read.csv("./Data/source_data/majors2008.1.csv")

# Read in geographic distance data
distdy <- read.csv("./Data/source_data/capdist.csv")

# Read in trade
strade <- read.csv("./Data/source_data/national_trade.csv")
dtrade <- read.csv("./Data/source_data/dyadic_trade.csv")

###### Manipulation for TERGM #######
## Times
# 1816 - 2000, without trade
# 1955 - 2000, with trade

# All alliance types
ally <- list()
allied <- atopdy$defense+atopdy$offense+atopdy$consul +atopdy$neutral
atop1 <- subset(atopdy,allied >0)
for(t in 1816:2000){
	cct <- subset(ssys, ssys$Year ==t)$CCode
	polyyr <- polity$year
	polityt <- subset(polity, polyyr ==t)
	cct <- polityt$ccode
	netmat <- matrix(0,length(cct),length(cct))
	atopt <- subset(atop1, atop1$year == t)
	cc1 <- atopt$mem1
	cc2 <- atopt$mem2
	ind1 <- match(cc1,cct)
	ind2 <- match(cc2,cct)
	netmat[cbind(ind1,ind2)] <- 1
	netmat[cbind(ind2,ind1)] <- 1
	ally[[t-1815]] <- netmat
}

# Just defense pacts
allyd <- list()
allied <- atopdy$defense
atop1 <- subset(atopdy,allied >0)
for(t in 1816:2000){
	cct <- subset(ssys, ssys$Year ==t)$CCode
	polyyr <- polity$year
	polityt <- subset(polity, polyyr ==t)
	cct <- polityt$ccode
	netmat <- matrix(0,length(cct),length(cct))
	atopt <- subset(atop1, atop1$year == t)
	cc1 <- atopt$mem1
	cc2 <- atopt$mem2
	ind1 <- match(cc1,cct)
	ind2 <- match(cc2,cct)
	netmat[cbind(ind1,ind2)] <- 1
	netmat[cbind(ind2,ind1)] <- 1
	allyd[[t-1815]] <- netmat
}

# Conflict History
library(combinat)
midhist <- list()
mids <- subset(maoz, maoz$dichmid20 ==1)
for(t in 1816:2000){
	cct <- subset(ssys, ssys$Year ==t)$CCode
	polyyr <- polity$year
	polityt <- subset(polity, polyyr ==t)
	cct <- polityt$ccode
	netmat <- matrix(0,length(cct),length(cct))
	midstm10 <- subset(mids,is.element(mids$year,(t-10):(t-1)))
	midscc1 <- paste(midstm10$statea,"_",midstm10$stateb,sep="")
	midscc2 <- paste(midstm10$stateb,"_",midstm10$statea,sep="")
	midscc <- c(midscc1,midscc2)
	alldy <- t(combn(cct,2))
	dycode <- paste(alldy[,1],"_",alldy[,2],sep="")
	dyccMid <- rbind(alldy[which(is.element(dycode,midscc)),])
	if(length(dyccMid)>0){
	ind1 <- match(dyccMid[,1],cct)
	ind2 <- match(dyccMid[,2],cct)
	netmat[cbind(ind1,ind2)] <- 1
	netmat[cbind(ind2,ind1)] <- 1
	}
	midhist[[t-1815]] <- netmat
}

# Common Enemy
comen <- list()
for(t in 1816:2000){
	cct <- subset(ssys, ssys$Year ==t)$CCode
	polyyr <- polity$year
	polityt <- subset(polity, polyyr ==t)
	cct <- polityt$ccode
	netmat <- matrix(0,length(cct),length(cct))
	midht <- midhist[[t-1815]]
	for(j in 2:length(cct)){
		for(h in 1:j){
			rj <- midht[j,-c(j,h)]
			rh <- midht[h,-c(j,h)]
			cejh <- as.numeric(sum(rj*rh)>0)
			netmat[j,h] <- cejh
			netmat[h,j] <- cejh
		}
	}			
	comen[[t-1815]] <- netmat
}


# Monadic Conflict History
library(combinat)
threat <- list()
mids <- subset(maoz, maoz$dichmid20 ==1)
for(t in 1816:2000){
	cct <- subset(ssys, ssys$Year ==t)$CCode
	polyyr <- polity$year
	polityt <- subset(polity, polyyr ==t)
	cct <- polityt$ccode
	netmat <- matrix(0,length(cct),length(cct))
	midstm10 <- subset(mids,is.element(mids$year,(t-10):(t-1)))
	midscc1 <- midstm10$statea
	midscc2 <- midstm10$stateb
	midscc <- c(midscc1,midscc2)
	for(i in 2:length(cct)){
			midsi <- length(which(midscc==cct[i]))
		for(j in 1:i){
				midsj <- length(which(midscc==cct[j]))
				netmat[i,j] <- midsi + midsj
				netmat[j,i] <- midsi + midsj
		}
	}
	threat[[t-1815]] <- netmat
}



trade <- list()
for(t in 1955:2000){
	cct <- subset(ssys, ssys$Year ==t)$CCode
	polyyr <- polity$year
	polityt <- subset(polity, polyyr ==t)
	cct <- polityt$ccode
	netmat <- matrix(0,length(cct),length(cct))
	dtradet <- subset(dtrade,dtrade$year==t)
	ind1 <- match(dtradet$ccode1,cct)
	ind2 <- match(dtradet$ccode2,cct)
	inds <- na.omit(cbind(ind1,ind2,dtradet$flow1,dtradet$flow2))
	netmat[inds[,c(1,2)]] <- inds[,3] + inds[,4]
	netmat[inds[,c(2,1)]] <- inds[,3] + inds[,4]		
	trade[[t-1815]] <- netmat
}


# Contiguity
contiguous <- list()
for(t in 1816:2000){
	cct <- subset(ssys, ssys$Year ==t)$CCode
	polyyr <- polity$year
	polityt <- subset(polity, polyyr ==t)
	cct <- polityt$ccode
	netmat <- matrix(0,length(cct),length(cct))
	ind1 <- match(contig$statelno,cct)
	ind2 <- match(contig$statehno,cct)
	ind1c <- ind1[which(complete.cases(cbind(ind1,ind2)))]
	ind2c <- ind2[which(complete.cases(cbind(ind1,ind2)))]
	netmat[cbind(ind1c,ind2c)] <- 1
	netmat[cbind(ind2c,ind1c)] <- 1
	contiguous[[t-1815]] <- netmat
}


# Joint Democracy
# Abs difference
jntdem <- list()
demdif <- list()
for(t in 1816:2000){
	cct <- subset(ssys, ssys$Year ==t)$CCode
	polyyr <- polity$year
	polityt <- subset(polity, polyyr ==t)
	cct <- polityt$ccode
	netmat1 <- matrix(0,length(cct),length(cct))
	netmat2 <- matrix(0,length(cct),length(cct))
	polyyr <- polity$year
	polityt <- subset(polity, polyyr ==t)
	for(j in 2:length(cct)){
		for(h in 1:j){
			demj <- polityt$democ[which(polityt$ccode==cct[j])]
			demh <- polityt$democ[which(polityt$ccode==cct[h])]
			if(length(demj)*length(demh)>0){
			netmat1[j,h] <- as.numeric(demj >= 6)*as.numeric(demh >= 6)
			netmat1[h,j] <- as.numeric(demj >= 6)*as.numeric(demh >= 6)
			netmat2[j,h] <- abs(demj-demh)
			netmat2[h,j] <- abs(demj-demh)
			}
		}
	}			
	jntdem[[t-1815]] <- netmat1
	demdif[[t-1815]] <- netmat2
}


# Column ccode
# rowccode
rowcc <- list()
colcc <- list()
for(t in 1816:2000){
	cct <- subset(ssys, ssys$Year ==t)$CCode
	polyyr <- polity$year
	polityt <- subset(polity, polyyr ==t)
	cct <- polityt$ccode
	netmatc <- matrix(0,length(cct),length(cct))
	netmatr <- matrix(0,length(cct),length(cct))
	for(i in 1:length(cct)){
		netmatc[,i] <- cct[i]
		netmatr[i,] <- cct[i]
	}
	rowcc[[t-1815]] <- netmatr
	colcc[[t-1815]] <- netmatc
}

# major powers
majpow <- list()
for(t in 1816:2000){
	cct <- subset(ssys, ssys$Year ==t)$CCode
	polyyr <- polity$year
	polityt <- subset(polity, polyyr ==t)
	cct <- polityt$ccode
	majpt <- subset(majp,majp$StYear <= t & majp$EndYear >= t)
	netmat <- matrix(0,length(cct),length(cct))
	for(i in 2:length(cct)){
		mpi <- is.element(cct[i],majpt$CCode)
		for(j in 1:i){
			mpj <- is.element(cct[j],majpt$CCode)
			netmat[i,j] <- mpi + mpj
			netmat[j,i] <- mpi + mpj
		}
	}
	majpow[[t-1815]] <- netmat
}


# dyadic distance
dydist <- list()
for(t in 1816:2000){
	cct <- subset(ssys, ssys$Year ==t)$CCode
	polyyr <- polity$year
	polityt <- subset(polity, polyyr ==t)
	cct <- polityt$ccode
	ind1 <- match(distdy$numa,cct)
	ind2 <- match(distdy$numb,cct)
	distdyt <- na.omit(cbind(ind1,ind2,distdy))
	netmat <- matrix(0,length(cct),length(cct))
	netmat[,as.matrix(distdyt[,c(1,2)])] <- distdyt$midist
	netmat[,as.matrix(distdyt[,c(2,1)])] <- distdyt$midist
	dydist[[t-1815]] <- netmat
}


#### Make datasets ####

## all alliances, 1816 ##

all1816 <- NULL
for(t in 1816:2000){
	i <- t - 1815
	allyt <- ally[[i]]
	midhistt <- midhist[[i]]
	coment <- comen[[i]]
	threatt <- threat[[i]]
	contigt <- contiguous[[i]]
	jntdemt <- jntdem[[i]]
	demdift <- demdif[[i]]
	majpowt <- majpow[[i]]
	dydistt <- dydist[[i]]
	rowcct <- rowcc[[i]]
	colcct <- colcc[[i]]
	mpledat <- ergmMPLE(network(allyt,directed=F)~edges+kstar(2)+gwesp(1.25)+edgecov(midhistt)+edgecov(coment)+edgecov(threatt)+edgecov(contigt)+edgecov(jntdemt)+edgecov(demdift)+edgecov(majpowt)+edgecov(majpowt)+edgecov(dydistt)+edgecov(rowcct)+edgecov(colcct), fitmodel=F)
	datt <- cbind(mpledat$response,mpledat$predictor,t)
	all1816 <- rbind(all1816,datt)
}
all1816 <- data.frame(all1816)
nm <- names(all1816)
spnm <- strsplit(nm,split="\\.")
nam <- NULL
for(i in 1:length(spnm)){
	nam <- c(nam,spnm[[i]][[min(c(2,length(spnm[[i]])))]])
}
nam[1] <- "y"
names(all1816) <- nam



## all alliances, 1955##
all1955 <- NULL
for(t in 1955:2000){
	i <- t - 1815
	allyt <- ally[[i]]
	midhistt <- midhist[[i]]
	coment <- comen[[i]]
	threatt <- threat[[i]]
	contigt <- contiguous[[i]]
	jntdemt <- jntdem[[i]]
	demdift <- demdif[[i]]
	tradet <- trade[[i]]
	majpowt <- majpow[[i]]
	dydistt <- dydist[[i]]
	rowcct <- rowcc[[i]]
	colcct <- colcc[[i]]
	mpledat <- ergmMPLE(network(allyt,directed=F)~edges+kstar(2)+gwesp(1.25)+edgecov(midhistt)+edgecov(coment)+edgecov(threatt)+edgecov(contigt)+edgecov(jntdemt)+edgecov(demdift)+edgecov(majpowt)+edgecov(majpowt)+edgecov(dydistt)+edgecov(tradet)+edgecov(rowcct)+edgecov(colcct), fitmodel=F)
	datt <- cbind(mpledat$response,mpledat$predictor,t)
	all1955 <- rbind(all1955,datt)
}
all1955 <- data.frame(all1955)
nm <- names(all1955)
spnm <- strsplit(nm,split="\\.")
nam <- NULL
for(i in 1:length(spnm)){
	nam <- c(nam,spnm[[i]][[min(c(2,length(spnm[[i]])))]])
}
nam[1] <- "y"
names(all1955) <- nam

## def alliances, 1816 ##

def1816 <- NULL
for(t in 1816:2000){
	i <- t - 1815
	allyt <- allyd[[i]]
	midhistt <- midhist[[i]]
	coment <- comen[[i]]
	threatt <- threat[[i]]
	contigt <- contiguous[[i]]
	jntdemt <- jntdem[[i]]
	demdift <- demdif[[i]]
	majpowt <- majpow[[i]]
	dydistt <- dydist[[i]]
	rowcct <- rowcc[[i]]
	colcct <- colcc[[i]]
	mpledat <- ergmMPLE(network(allyt,directed=F)~edges+kstar(2)+gwesp(1.25)+edgecov(midhistt)+edgecov(coment)+edgecov(threatt)+edgecov(contigt)+edgecov(jntdemt)+edgecov(demdift)+edgecov(majpowt)+edgecov(majpowt)+edgecov(dydistt)+edgecov(rowcct)+edgecov(colcct), fitmodel=F)
	datt <- cbind(mpledat$response,mpledat$predictor,t)
	def1816 <- rbind(def1816,datt)
}
def1816 <- data.frame(def1816)
nm <- names(def1816)
spnm <- strsplit(nm,split="\\.")
nam <- NULL
for(i in 1:length(spnm)){
	nam <- c(nam,spnm[[i]][[min(c(2,length(spnm[[i]])))]])
}
nam[1] <- "y"
names(def1816) <- nam



## all alliances, 1955##
def1955 <- NULL
for(t in 1955:2000){
	i <- t - 1815
	allyt <- allyd[[i]]
	midhistt <- midhist[[i]]
	coment <- comen[[i]]
	threatt <- threat[[i]]
	contigt <- contiguous[[i]]
	jntdemt <- jntdem[[i]]
	demdift <- demdif[[i]]
	tradet <- trade[[i]]
	majpowt <- majpow[[i]]
	dydistt <- dydist[[i]]
	rowcct <- rowcc[[i]]
	colcct <- colcc[[i]]
	mpledat <- ergmMPLE(network(allyt,directed=F)~edges+kstar(2)+gwesp(1.25)+edgecov(midhistt)+edgecov(coment)+edgecov(threatt)+edgecov(contigt)+edgecov(jntdemt)+edgecov(demdift)+edgecov(majpowt)+edgecov(majpowt)+edgecov(dydistt)+edgecov(tradet)+edgecov(rowcct)+edgecov(colcct), fitmodel=F)
	datt <- cbind(mpledat$response,mpledat$predictor,t)
	def1955 <- rbind(def1955,datt)
}
def1955 <- data.frame(def1955)
nm <- names(def1955)
spnm <- strsplit(nm,split="\\.")
nam <- NULL
for(i in 1:length(spnm)){
	nam <- c(nam,spnm[[i]][[min(c(2,length(spnm[[i]])))]])
}
nam[1] <- "y"
names(def1955) <- nam


write.csv(all1816,"./Data/all1816.csv")
write.csv(all1955,"./Data/all1955.csv")
write.csv(def1816,"./Data/def1816.csv")
write.csv(def1955,"./Data/def1955.csv")

### Run the Models to get MPLE estimates ###
modFC <- glm(y ~ kstar2 + gwesp + midhistt + coment + threatt + contigt + jntdemt + demdift + majpowt + dydistt,family=binomial, data=all1816)
modFN <- glm(y ~ kstar2 + gwesp + midhistt + coment + threatt + contigt + jntdemt + demdift + majpowt + dydistt,family=binomial, data=def1816)
modHC <- glm(y ~ kstar2 + gwesp + midhistt + coment + threatt + contigt + jntdemt + demdift + majpowt + dydistt+tradet,family=binomial, data=all1955)
modHN <- glm(y ~ kstar2 + gwesp + midhistt + coment + threatt + contigt + jntdemt + demdift + majpowt + dydistt+tradet,family=binomial, data=def1955)

lmodFC <- glm(y ~ midhistt + coment + threatt + contigt + jntdemt + demdift + majpowt + dydistt,family=binomial, data=all1816)
lmodFN <- glm(y ~ midhistt + coment + threatt + contigt + jntdemt + demdift + majpowt + dydistt,family=binomial, data=def1816)
lmodHC <- glm(y ~ midhistt + coment + threatt + contigt + jntdemt + demdift + majpowt + dydistt+tradet,family=binomial, data=all1955)
lmodHN <- glm(y ~ midhistt + coment + threatt + contigt + jntdemt + demdift + majpowt + dydistt+tradet,family=binomial, data=def1955)

### Make ROC Plots ####
predFC <- prediction(predict(modFC, type="resp"),modFC$y)
aucFC <- attr(performance(predFC,measure="auc"),"y.values")[[1]]
rocFC <- performance(predFC,measure="tpr",x.measure = "fpr")

predFN <- prediction(predict(modFN, type="resp"),modFN$y)
aucFN <- attr(performance(predFN,measure="auc"),"y.values")[[1]]
rocFN <- performance(predFN,measure="tpr",x.measure = "fpr")

predHC <- prediction(predict(modHC, type="resp"),modHC$y)
aucHC <- attr(performance(predHC,measure="auc"),"y.values")[[1]]
rocHC <- performance(predHC,measure="tpr",x.measure = "fpr")

predHN <- prediction(predict(modHN, type="resp"),modHN$y)
aucHN <- attr(performance(predHN,measure="auc"),"y.values")[[1]]
rocHN <- performance(predHN,measure="tpr",x.measure = "fpr")

predlFC <- prediction(predict(lmodFC, type="resp"),lmodFC$y)
auclFC <- attr(performance(predlFC,measure="auc"),"y.values")[[1]]
roclFC <- performance(predlFC,measure="tpr",x.measure = "fpr")

predlFN <- prediction(predict(lmodFN, type="resp"),lmodFN$y)
auclFN <- attr(performance(predlFN,measure="auc"),"y.values")[[1]]
roclFN <- performance(predlFN,measure="tpr",x.measure = "fpr")

predlHC <- prediction(predict(lmodHC, type="resp"),lmodHC$y)
auclHC <- attr(performance(predlHC,measure="auc"),"y.values")[[1]]
roclHC <- performance(predlHC,measure="tpr",x.measure = "fpr")

predlHN <- prediction(predict(lmodHN, type="resp"),lmodHN$y)
auclHN <- attr(performance(predlHN,measure="auc"),"y.values")[[1]]
roclHN <- performance(predlHN,measure="tpr",x.measure = "fpr")

### Keep only model summary for space ###
modFC <- summary(modFC)
modFN <- summary(modFN)
modHC <- summary(modHC)
modHN <- summary(modHN)
lmodFC <- summary(lmodFC)
lmodFN <- summary(lmodFN)
lmodHC <- summary(lmodHC)
lmodHN <- summary(lmodHN)
gc()

### ROC Plots ###
pdf("roc_all_full.pdf",height=6, width=6.5, pointsize=12, family="Times")
par(mar=c(4,4,1,1),las=1, cex.lab=1.5, col.axis="white")
plot(rocFC, xaxt="n", yaxt="n",xlab = "% Non-Alliances Predicted to be Alliances", ylab="% Alliances Predicted to be Alliances", lwd=2.5)
plot(roclFC, add=T, lwd=2.5, col="grey65", xaxt="n", yaxt="n")
axis(1, at=seq(0,1,by=.2), label = 100*seq(0,1,by=.2), col.axis="black")
axis(2, at=seq(0,1,by=.2), label = 100*seq(0,1,by=.2), col.axis="black")
legend("bottomright", legend=c(paste("TERGM",", AUC = ",round(aucFC,3),sep=""),paste("Logit",", AUC = ",round(auclFC,3),sep="")), lwd=2.5, col=c("black", "grey65"), cex=1.5)
dev.off() 

pdf("roc_def_full.pdf",height=6, width=6.5, pointsize=12, family="Times")
par(mar=c(4,4,1,1),las=1, cex.lab=1.5, col.axis="white")
plot(rocFN, xaxt="n", yaxt="n",xlab = "% Non-Alliances Predicted to be Alliances", ylab="% Alliances Predicted to be Alliances", lwd=2.5)
plot(roclFN, add=T, lwd=2.5, col="grey65", xaxt="n", yaxt="n")
axis(1, at=seq(0,1,by=.2), label = 100*seq(0,1,by=.2), col.axis="black")
axis(2, at=seq(0,1,by=.2), label = 100*seq(0,1,by=.2), col.axis="black")
legend("bottomright", legend=c(paste("TERGM",", AUC = ",round(aucFN,3),sep=""),paste("Logit",", AUC = ",round(auclFN,3),sep="")), lwd=2.5, col=c("black", "grey65"), cex=1.5)
dev.off() 

pdf("roc_all_half.pdf",height=6, width=6.5, pointsize=12, family="Times")
par(mar=c(4,4,1,1),las=1, cex.lab=1.5, col.axis="white")
plot(rocHC, xaxt="n", yaxt="n",xlab = "% Non-Alliances Predicted to be Alliances", ylab="% Alliances Predicted to be Alliances", lwd=2.5)
plot(roclHC, add=T, lwd=2.5, col="grey65", xaxt="n", yaxt="n")
axis(1, at=seq(0,1,by=.2), label = 100*seq(0,1,by=.2), col.axis="black")
axis(2, at=seq(0,1,by=.2), label = 100*seq(0,1,by=.2), col.axis="black")
legend("bottomright", legend=c(paste("TERGM",", AUC = ",round(aucHC,3),sep=""),paste("Logit",", AUC = ",round(auclHC,3),sep="")), lwd=2.5, col=c("black", "grey65"), cex=1.5)
dev.off() 

pdf("roc_def_half.pdf",height=6, width=6.5, pointsize=12, family="Times")
par(mar=c(4,4,1,1),las=1, cex.lab=1.5, col.axis="white")
plot(rocHN, xaxt="n", yaxt="n",xlab = "% Non-Alliances Predicted to be Alliances", ylab="% Alliances Predicted to be Alliances", lwd=2.5)
plot(roclHN, add=T, lwd=2.5, col="grey65", xaxt="n", yaxt="n")
axis(1, at=seq(0,1,by=.2), label = 100*seq(0,1,by=.2), col.axis="black")
axis(2, at=seq(0,1,by=.2), label = 100*seq(0,1,by=.2), col.axis="black")
legend("bottomright", legend=c(paste("TERGM",", AUC = ",round(aucHN,3),sep=""),paste("Logit",", AUC = ",round(auclHN,3),sep="")), lwd=2.5, col=c("black", "grey65"), cex=1.5)
dev.off() 

## Remove data for memory ##

rm(all1816)
rm(def1816)
rm(all1955)
rm(def1955)
gc()

### Run bootstrap TERGM ###
source("./Code/bootstrap_tergm.R")

### Read in bootstrapped results ####
resall1816 <- NULL
for(i in 1:5){
	filei <- paste("./Data/all1816res",i,".csv",sep="")
	resall1816 <- rbind(resall1816,read.csv(filei))
}

resdef1816 <- NULL
for(i in 1:5){
	filei <- paste("./Data/def1816res",i,".csv",sep="")
	resdef1816 <- rbind(resdef1816,read.csv(filei))
}

resall1955 <- NULL
for(i in 1:5){
	filei <- paste("./Data/all1955res",i,".csv",sep="")
	resall1955 <- rbind(resall1955,read.csv(filei))
}

resdef1955<- NULL
for(i in 1:5){
	filei <- paste("./Data/def1955res",i,".csv",sep="")
	resdef1955 <- rbind(resdef1955,read.csv(filei))
}


### Make Results Tables ###
## First all alliance types
rown <- rownames(modHC$coef)

estFC <- rep(NA,length(rown))
seFC <- rep(NA,length(rown))
estFC[match(rownames(modFC$coef),rown)] <- modFC$coef[,1]
resall1816 <- resall1816[,2:ncol(resall1816)]
names(resall1816)[1] <- rown[1]
seFC[match(names(resall1816),rown)] <- apply(resall1816,2,sd)

estHC <- rep(NA,length(rown))
seHC <- rep(NA,length(rown))
estHC[match(rownames(modHC$coef),rown)] <- modHC$coef[,1]
resall1955 <- resall1955[,2:ncol(resall1955)]
names(resall1955)[1] <- rown[1]
seHC[match(names(resall1955),rown)] <- apply(resall1955,2,sd)

estlFC <- rep(NA,length(rown))
selFC <- rep(NA,length(rown))
estlFC[match(rownames(lmodFC$coef),rown)] <- lmodFC$coef[,1]
selFC[match(rownames(lmodFC$coef),rown)] <- lmodFC$coef[,2]

estlHC <- rep(NA,length(rown))
selHC <- rep(NA,length(rown))
estlHC[match(rownames(lmodHC$coef),rown)] <- lmodHC$coef[,1]
selHC[match(rownames(lmodHC$coef),rown)] <- lmodHC$coef[,2]

resmatall <- cbind(estlFC,selFC,estFC,seFC,estlHC,selHC,estHC,seHC)
rownames(resmatall) <- rown


## Defense alliances only
rown <- rownames(modHC$coef)

estFN <- rep(NA,length(rown))
seFN <- rep(NA,length(rown))
estFN[match(rownames(modFN$coef),rown)] <- modFN$coef[,1]
resdef1816 <- resdef1816[,2:ncol(resdef1816)]
names(resdef1816)[1] <- rown[1]
seFN[match(names(resdef1816),rown)] <- apply(resdef1816,2,sd)

estHN <- rep(NA,length(rown))
seHN <- rep(NA,length(rown))
estHN[match(rownames(modHN$coef),rown)] <- modHN$coef[,1]
resdef1955 <- resdef1955[,2:ncol(resdef1955)]
names(resdef1955)[1] <- rown[1]
seHN[match(names(resdef1955),rown)] <- apply(resdef1955,2,sd)

estlFN <- rep(NA,length(rown))
selFN <- rep(NA,length(rown))
estlFN[match(rownames(lmodFN$coef),rown)] <- lmodFN$coef[,1]
selFN[match(rownames(lmodFN$coef),rown)] <- lmodFN$coef[,2]

estlHN <- rep(NA,length(rown))
selHN <- rep(NA,length(rown))
estlHN[match(rownames(lmodHN$coef),rown)] <- lmodHN$coef[,1]
selHN[match(rownames(lmodHN$coef),rown)] <- lmodHN$coef[,2]

resmatdef <- cbind(estlFN,selFN,estFN,seFN,estlHN,selHN,estHN,seHN)
rownames(resmatdef) <- rown

library(xtable)
xtable(resmatall,digit=7)
xtable(resmatdef,digit=7)




